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RECONSTRUCTION OF LAME MODULI AND DENSITY AT THE 
BOUNDARY ENABLING DIRECTIONAL ELASTIC WAVEFIELD 

DECOMPOSITION 

MAARTEN V. DE HOOP*, GEN NAKAMURAt, AND JIAN ZHAI* 

Abstract. We consider the inverse boundary value problem for the system of equations describ¬ 
ing elastic waves in isotropic media on a bounded domain in via a finite-time Laplace transform. 
The data is the dynamical Dirichlet-to-Neumann map. More precisely, using the full symbol of the 
transformed Dirichlet-to-Neumann map viewed as a semiclassical pseudodifferential operator, we give 
an explicit reconstruction of both Lame parameters and the density, as well as their derivatives, at 
the boundary. We also show how this boundary reconstruction leads to a decomposition of incoming 
and outgoing waves. 
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1. Introduction. 

We let n C be a bounded domain with a smooth boundary 50. We consider 
the following initial boundary value problem for the system of equations describing 
elastic waves 

{ pd^u = div(Ce(M)) =: Lu in VLt = O x (0,T), 
u=f onE = 50x (0,T), 
u{x, 0) = dtu{x, 0) = 0 in O, 

with f{x,0) = 0 and ^/(x, 0) = 0 for x € dfl. Here, u denotes the displacement 
vector and e{u) = {eij{u)) = (Vm + (Vu)^)/2 the linear strain tensor which is the 
symmetric part of Vm. Furthermore, C = C(x) = (Cijki(x)) is the elasticity tensor 
and p is the density of mass. We assume that C is isotropic, that is, 

(1.2) Cijkl^x') = \{x^5ij 5kl “t“ P'i^X^i^SikSjl -f SnSjk) 

with Kronecker’s delta Sij and Lame moduli X, p, G such that p > 0 and 

A + 2/i > 0 on O. Also, p G and p > 0 on O. 

The hyperbolic or dynamical Dirichlet-to-Neumann map (DN map) At is defined 
according to 

(1.3) At : R2(E) 9 / ^ Olu := (Ce(M))i^|an £ U([0, T], 

where u is the solution of (1.1), Ce(u) is a 3 x 3 matrix with its (*, fc) component 
{Ce{u))ik given by {Ce{u))ik = v is the outward unit normal to 

dfl. Physically, Olu denotes the traction at dfl. 

In this paper, we consider the inverse problem of recovering A, p, p, as well as 
all their derivatives, at the boundary 9D from At. Our major result for this inverse 
problem is as follows. 
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Reconstruction of Lame parameters and density at the boundary 


Theorem 1.1. The DN map Ax identifies X, p and all their derivatives on dTl 
uniquely. There is an explicit reconstruction procedure for these identification. 

Remark 1.2. Since the procedure we present to recover X, p, p and their deriva¬ 
tives at dtd is local, we also have a localized version of Theorem 1.1 with partial bound¬ 
ary data. That is, in the definition of Ax we can replace (Ce(u))j^| 9 Q by (Ce(w))t^|ro 
and confine the boundary sources, f, to those with supp f{-,l) C Fp (t G (0,r)), where 
Fp is a relatively open subset of dVL. We can recover X,p,p and all their derivatives 
ai Fp. As an additional explanation for this which should be given later, see the last 
paragraph of Section 2. 

The uniqueness of the inverse problem considered here was established by Rachele 
[16], but giving a procedure to recover the parameters, (A,/i, p), has been left open 
for over 15 years. One of the complications is the occurrence of two metrics in the 
dynamical system of equations that cannot be straightforwardly separated at the 
boundary. Indeed, the usual special solutions including high-frequency asymptotic 
ones or progressive wave solutions based on polarization decoupling are coupled at 
the boundary. 

Moreover, the determination of (A, p, p) on boundary implies the following 

Corollary 1.3. In addition to the conditions appearing in Theorem 1.1, let 
X, p, p be real analytic in the neighborhood of II. Then Ax determines uniquely X,p,p 
on n. 

A brief remark for this corollary should be given. Although there is a more general 
result by Rachele [16],[17], the context and argument of deriving this corollary differs 
from those of Rachele’s. 

For the static elastic inverse boundary value problem, an explicit reconstruction 
of A and p at the boundary from the full symbol of the static DN map was obtained 
[13, 14]. For the reconstruction of a transversely isotropic elasticity tensor, see [15]. 
The approach was originally developed by Sylvester and Uhlmann [21] for the elec¬ 
trical impedance tomography problem. The approach is also applied to Maxwell’s 
equations [20, 12]. We generalize this type of reconstruction to dynamical elastic 
inverse boundary value problems. We note that our procedure is quite general and 
can be extended from isotropy to anisotropy with certain symmetries, which is the 
subject of a forthcoming paper. 

The key component of the reconstruction is a connection between Ax and the 
asymptotic expansion of the DN map, a semiclassical pseudodifferential operator A^ 
say, for some elliptic system of equations containing a small parameter h via a finite¬ 
time Laplace transform. M. Ikehata has been using the finite Laplace transform ef¬ 
fectively to develop his enclosure method both for parabolic equations and hyperbolic 
equations, see [8] and reference therein. For the convenience of our description, the 
partial differential operator of this system is referred by M. We will identify (A, p, p) 
and all of their derivatives from A^ by factorizing Ai into the product of two first 
order semiclassical pseudodifferential operators with small parameter /i = i, where r 
is nothing but the Laplace variable of this transform. Also this factorization is noth¬ 
ing but the one used to provide the up/down going decomposition of waves which is 
equivalent to the incoming/downgoing decomposition of waves in the Laplace domain. 
Further this decomposition can be linked to the corresponding decomposition in the 
space time domain. We will briefly discuss about this connection more precisely at 
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the end of this paper. The up/down going decomposition for scalar waves is discussed 
in [19], and for elastic waves in [6]. In this paper, we connect the up/down going 
decomposition with Dirichlet-to-Neumann map. 



Fig. 1. incoming/outgoing waves. 

Concerning our inverse problem, there are two byproducts of the factorization. 
The one is the explicit form of the principal symbol of A^(s) and relation of its s- 
derivatives to the non-principal symbols of A^(s). Here A^ is the DN map defined 
likewise A^ on the boundary r(s) of the subdomain H(s) = {x € : dist(a:,90) > s} 

of O with 0 < s 1. The other is that a modification A(s) (cf. (3.20)) of A^(s) 
satisfies a Riccati type equation. By solving the Riccati type equation, which is an 
initial value problem, one can propagate the data into at least a thin layer near the 
boundary. This technique is known as invariant embedding in extensive geophysics 
literature [2, 4, 5, 6, 11]. 

Knowing all the derivatives of the coefficients at and using the Riccati equation 
for A(s), we can generate an approximation of At{s) on r(s). Then, we can get an 
approximation of A, /i, p and their derivatives at r(s). Repeating this process, leads to 
an approximation for A, /i, p layer by layer in the interior of H, using the DN map At as 
the data. The associated algorithm is called layer stripping. The layer stripping was 
first developed for the electrical impedance tomography problem in [3, 18]. Nakamura, 
Tanuma and Uhlmann [15] developed such an algorithm for the static elastic inverse 
boundary value problem in the case of transverse isotropy. 

The key application of the problem we are considering is (reflection) seismology. 
In actual seismic acquisition, raw vibroseis data are modeled by the local Neumann- 
to-Dirichlet (ND) map: The boundary values are given by the normal traction un¬ 
derneath the base plate of a vibroseis and are zero (free surface) elsewhere, while the 
particle displacement (in fact, velocity) is measured by geophones located in a subset 
of the boundary (Earths surface) (see [1]). Although, the local dynamical Dirichlet- 
to-Neumann map and Neumann-to-Dirichlet map do not have the same information, 
the transformed ND map and transformed DN map are microlocally inverse to each 
other. Since we are only dealing with the symbol of transformed DN map, the results 
of this paper apply to the practical setting. 

The remainder of this paper is organized as follows. In Section 2, we give an 
asymptotic identity which connects At with A^ via the finite-time Laplace transform 
introducing variable f = ^. Based on this identity, we only have to find a recon¬ 
struction procedure using A^. In a similar way, DN maps A^(s) are defined on each 
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r(s). In Section 3, the full symbol of seiniclassical pseudodifferential operator A^(s) 
is analyzed using the factorization of A4. Also, as a byproduct of the factorization, 
we develop a layer stripping algorithm. Section 4 is devoted to giving a procedure and 
formulas for the reconstruction of (A, /i, p) and all their derivatives, from the explicit 
form of the principal symbol of A^(s) at each r(s) for 0 < s <C 1 in terms of the 
boundary normal coordinates associated with r(s). In the final section, Section 5, we 
will discuss about an another implication of the factorization. That is we give the 
aforementioned link between the outgoing/incoming decomposition of waves in the 
Laplace domain and that of in the space time domain. 

2. Reduction to an elliptic boundary value problem with a small pa¬ 
rameter. 

First we introduce a family of symbol classes for semiclassical pseudodifferen¬ 
tial operators. Let x (0, hg) —>■ be a function that is smooth in 

(x, depending on h G (0, ho] with a small ho > 0. We say that for m G R, A 

belongs to a symbol class 5(m), if for any a, /3 G Z" , there exists a constant > 0 
such that 


\D^DlA{x,i;h)\ < (x,e) gM^", (0,/io], 

where Z+ = N U {0}, {^) = \Jl + |^p. 

We say that A G S{m) is called a classical symbol if for any a G Z” and N = 
1, 2, • • •, there exist a constant Ca^iv > 0 such that 


Af-l 

|5“(A- ^ h^A,)\ < C^,Nh^{^r, (x,e) G 

3=0 

with each Aj G S(m) independent of h, and write 


A ^ mod 0{h°°S{m)). 

3=0 

Aq and are called the the principal symbol cr(A) and full symbol h{A) of 

A, respectively. In this paper we only consider classical symbols. We also denote for 
A, A' G S{m), 


A ^ A' mod 0{h^S{m)), 
if A - A' = h^B with B G S{m). 

For this family of symbol classes S{m) with m G R, we use the standard theory 
of semiclassical pseudodifferential operators (see [10],[23] ). Also we use the above 
notations and terminologies for symbol classes associated to semiclassical pseudod¬ 
ifferential operators on compact manifolds. In the sequel Q will be either g = 3 or 
q = 6 which can be easily noticed in contexts. Furthermore, we abuse the notation 
0{h^S{m)) to use it also for the associated semiclassical pseudodifferential operators. 

Now we consider the following boundary value problem: 


( 2 . 1 ) 


pv — div(Ce:(u)) =0 in fl, 
V = ip on dfl 
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with a small (real-valued) parameter h G (0,ho]. We define the corresponding DN 
map for (2.1) according to 

J^h . ^ ^ ^ h{Ce{v))iy\on £ 


where v solves (2.1). In a likewise fashion, we define A^(s) by replacing fl by n(s) 
and dfl by r(s) while replacing f hy i/j G iJ®/^(r(s)) (cf. (3.1)) emphasizing the s 
dependence. Of course, = A^(0). We note that A^ and A^(s) are semiclassical 
pseudodifferential operators belonging to the class 5(1) with n = 2. 

We show that we can obtain the full symbol of A^ from A^ via a finite-time 
Laplace transform. First we introduce the finite-time Laplace transform w G 

oiuGC{[0,T],H^mby 


w{x,t) = [Ctu){x,t) 


L[x,t)e "^‘dt with T > 0. 


Jo 

In order to establish the connection between and A, we let x(<) = 
and define 


t2(t e [o,r]) 


At : W((0,r);5O) 


by 


Krpcj) = Kt{x4')- 


For any (j) G let u solve (1.1) with boundary value / = x4>- By the 

estimates for solutions of hyperbolic system (1.1), we have 

\\9iu{-,T)\\fj2-^Q) = 0 {\\x4>\\h^{^)) = ^{\\4>\\m/'^{dn)) 
for j = 0,1 (see [22]). Becanse 

Crid^u) = r’^Criu) + dtu{T)e~^^ + ru(T)e“'^^ — dtu{0) — tu(0), 
we have that by applying £t to both sides of (1.1), and divide by = -^ 


( 2 . 2 ) 


pw — h^div(Ce(w)) = r in fl, 
w = C-t{x4') on dVL, 


where w = Ct(u) and r has an estimate ||r||^i(f 2 ) < C'e“'‘'^^||(/)||^5/2(9f2) with some 
constant C for any given k satisfying 0 < k < 1. 

Subtracting (2.1) from (2.2) with tp = Ct{x4’)i find that z = w — v satisfies 


(2.3) 

Hence, we have 


pz — div(Ce(z)) = r in H, 
u = 0 on on. 


Ct {Olu) = Olw = Olv + Blz, 


with 


11 2 II ^3/2 ( 90 ) < C||r||^^i(o) < Ce 11^5/2 (go), 
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by standard elliptic regularity theory. We observe that 

hdLV = hA'^(CT{x(l>)) 


and 


hCr {dhu) = hCTAxix'f)- 

Thus, defining Ct ■ —?> H^/'^{dVL) by Criij)) = ^^t{x4')i can rewrite the 

formula above as 


HC.t'Aj’ = £,x + 0{g 

Here, denotes an operator from to with the estimate 

||C>(e < Ce 

We note that Lt is just a multiplication by hence it is invertible and 

can be estimated by 0{t~^) for r 1. Therefore, 

modulo an operator in —>■ with the estimate 0{h°°). Therefore, 

we can obtain the full symbol of from CtAtC^^, due to what we will mention in 
the last paragraph of this section. 

We note that in the above one can choose any smooth function for x that is 
consistent with the initial conditions such that x(t)e“'^*dt behaves polynomially 
in T. 

The full symbol of a semiclassical pseudodifferential operator can be evaluated by 
applying to locally supported rapidly oscillating functions [23]. So the analysis can 
be local, and thus we have Remark 1.2. 

3. Analysis of the symbol of A^(s). 

Given a boundary point po € r(s), for any x € fl(s) near po, we use the boundary 
normal coordinates x = (x^(p),x^(p),x^) = (y^,y^,x^) = (y',x^), where p € r(s) is 
the nearest point to x, x^ = dist(x,p), and (x^ (p), x^ (p)) are the local coordinates 
of r(s) near po- Then r(s) is locally given as x^ = s. Let ('Ci)C 2 )C 3 )) be 

conormal vectors with respect to the coordinates (xi,X 2 ,X 3 ), such that 

- We will use the notation y = (r?i, 772 , 773 ) = {v',r] 3 )- In 
boundary normal coordinates, equation (2.1) attains the form 
(3.1) 

{ {MvY = pg^^Vk ^ = 0 in {x^ > s} for 1 < i < 3, 

v''\x3=s = 'Y\ 1 < 7 < 3. 

where Vj is the covariant derivative with respect to g|j and eki{v) = 2~^{yiVk+VkVi) 
is the linear strain tensor. 


Cb'='(x) 


E 


i9x* dx^ dx’^ dx^ ■ 

7 ^ 7^ ^abcd ) 

OXa uXb OXc OXd 
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with Cabcdix) given by (1.2). The induced metric G{x) = (g“(a;)) is given by 


r—1 


In terms of Jacobi matrix J = {dx°'/dxr; ^ < a, r < 3), G takes the form G = JJ^. 
The expression for A^(s) in boundary normal coordinates is 


(3.2) 


(A^(s)V')* = -^ X] 1 <i <3 


k.l=l 


at = s. The full symbol of A^(s) can be expanded as 

CT(A^(s))(y', rj') - X , rj') mod 0{h°^S{l)), 

j<0 

where each Xj{s){y' ,r]') G 5(1). 

Now, we define 


Q{x,v') = X 1 < b ^ < 3 , 


(3.3) 


2 


R{x,y') = 1 < b fc < 3j , 

D{x) = (C'*3'=3(x); l<i,k<3). 

The principal symbol, M{x,r]), of M is then given by 

(3.4) M{x, rf) = D{x)r]l + {R{x, rj') + R^{x, r]'))r ]3 + Q{x, rj') + p{x)G{x). 

By the assumption, M(x,r]) is a positive definite matrix for x G SI, 77 G 
Hence, for fixed (x, 77 '), det D~^Gm{x, = 0 in 773 admits 3 roots 773 = Q {j = 

1, 2,3) with positive imaginary parts and 3 roots Q (j = 1, 2, 3) with negative imagi¬ 
nary parts. Thus, 


Lemma 3.1 ([7]). There is a unique factorization 

M{x,r]) = H(x)“^/ 2 M(x, 77 )£)(x )"^/2 = (773 - 5^(x, 77'))(773 - 5 o(x, 77 ')), 

with Spec(5o(x, 77 ')) C C+, where Spec(5o(x, 77 ')) is the spectrum of So{x,ri'). In the 
above, 

So{x,r]') := CM{x,'n',()~^dC^ M(x, 77 ', C)"^dC 

where 7 C C+ := {^ G C : SC ■= imaginary part of C > 0} 7 S a continuous curve 
enclosing all the Q (j = 1,2,3). 

Then we have the following factorization of M{x,ri): 

M{x,ri) = (773 - 5 d(x, 77 '))L>(x )(773 - Soix,r]')), 



(3.5) 
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where 

So{x,rj') = D-^/\x)S^{x,ri')D^/\x). 

We arrive at 

Lemma 3.2. The operator M. admits a factorization 

. . M = {hDs - S*{x,hDy,;h) 

+hK{x,hDyi))D{x) {hDs — S{x,hDyi;h)), 

where S{x,ri'\h) G 5(1), K{x,'r]';h) G 5(0). Moreover, hDyi = {hDyi,hDy 2 ), hDyj 
—ihd/dy^ {j = 1,2) and the principal symbol, So{x,r]'), of S satisfies 

(3.7) Spec(5o(a;,?7')) c C+. 

Proof. Following (3.4), we write the full symbol a(A4) of Ai in the form, 

(3.8) a(M) = D{x)ril + [R{x, y) + {x, 77 '))%+ 

Q(x,r]') + p{x)G{x) + hFo{x)ri 3 + hFi{x,ij'), 
where Fi{x,ri') G 5(1), and Fq(x) is a matrix multiplication. We expand 

(hDs — S*{x, hDyi; h) + hK{x, hDyi\ h))D{x){hDs — S{x, hDy>; h)), 

yielding 

(3.9) h{DsD){hDs) - S*D{hDs) + hKD{hDs) - h{D,D)S + S*DS 

- hKDS - hD{DsS) - DS{hDs) + D{hDsf. 
Comparing (3.8) and (3.9), we find that 5 and K should satisfy 

(3.10) 

— S*{x, hDyi; h)D{x) + hK{x, hDyi; h)D{x) — D{x)S{x, hDyi; h) + hDsD{x) 

= R{x, hDyi) + R^{x, hDyi) + hFo{x) 

and 

(3.11) — h{DsD{x))S{x, hDyi] h) + S*(x, hDyi] h)D{x)S{x, hDyi] h) 

— hK{x, hDyi] h)D{x)S{x, hDyi] h) — hD{x){DsS{x, hDyi] h)) 

= hFi{x, hDyi) + Q{x, hDyi) + pG{x). 
Eliminating K in (3.11) by using (3.10), we get 

(3.12) ihDs)S + 52 + D-\R + R^ + hFo)S + hD-^Fi + D-^Q + D-^pG = 0. 
By the composition formula for symbols of pseudodifferential operators, we have 

(3-13) E p'] h)D!^,S{x, vf] h) 

a>0 

+ ^-^h^°‘^D^'{D~^{x){R{x,p') + R^{x,p'))D^,S{x,p']h) 

a>0 

+ hD~^{x)Fo{x)S{x,p']h) + hD~^{x)Fi{x,p') + D~^{x)Q{x,p') 


+ D ^{x)pG + hDsS{x,p']h) = 0. 
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We introduce the expansions 

S{x, 77 '; h) h~^Sj(x, T]'), 

3<0 

K{x,r)';h) Kj{x,t]') 

j<0 

with Sj € S(j + 1), Kj € S{j) for every j. We construct S via arranging terms of the 
same degree of h in (3.13). The terms of order 0(hP) give 

(3.14) D~^{x){R{x,r]') + {x,ri'))So{x,r]') + D~^{x)Q{x,ri') 

+ D~^{x)p{x)G + S'o(x, 77 ') = 0. 

Indeed, So (cf. (3.5)) satishes this equation. The terms of order 0{h) give 

(3.15) + D-\R + R^)S-i + ^D^'SoD^,So 

\a\ = l 

+ Y, ^D;,iD-\R + R^))DpSo + D-^FoSo + D-^F^ + D^Sq = 0 . 

I«l=i 

The terms which are of homogeneity of order 0{h~^) for j < —2 yield 

il'^l 

(3.16) SoSj+SjSo + D-^iR + R^)Sj+ Y 

l + m=j + \a\ 

\c.\>l 

+ Y + I] iD^'iD-\R + R^))D^,S,+i + D-^FoS,+i + D,Sj+i = 0. 

l“l=i 

To confirm that (3.15) and (3.16) have solutions, we note that 

SjSo + SoSj + D-\R + R^)Sj = -D-\Q + pG)S^^Sj + 
using (3.14). Since 

Spec(S'o) C C+, Spec(-£|-^(Q + pG)Sq^) C C+, 

we can indeed solve for Sj {j < —1) in (3.15) and (3.16). 

After constructing the full symbol of S, we determine the full symbol of K from 
(3.10). □ 

Proposition 3.3. Let XQ{s){y,p') be the principal symbol o/A^(s). Then 

(3.17) Ao(s)( 2 /', 77 ') = -i{D{x)So{x, 77 ') + R^{x, r]'))\^3^^. 


Proof. For a given s (0 < s <C 1), let T be given such that 0 < T — s -C 1. The 
parametrix U = U{y',x^;h) to the boundary value problem (2.1) satisfies locally 

MU - 0 mod 0{h°^S{2)) in x [s,7l, 

U\^3^j- = L 
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Equation (3.7) implies that the solution operator of the factor {hDs)—S*{x, hDyi; h)+ 
hK{x, hDy'] h) in the factorization (3.6), for decreasing s, is decaying of order 0{h°°). 
Hence, U satisfies 


{{hDs) - S{x, hDy.-h))U - 0 mod 0{h°^S{l)). 

Thus, 

{hDs)U\^^=s ^ S{x,hDyi\h)U\^3=s mod C>(h°°5(l)). 

Therefore, 

(3.18) 

A^{s){y', hDyf, h) ^ —i{D(x)S{x, hDy>; h) + (x, hDy/))\^3^g mod C>(/i°°5(l)). 

Then formula (3.17) follows immediately. □ 

Next we establish the relation between the principal symbol Ao(s) and the lower 
order ones Aj(s), j < —1. Below, we use the notation mod {Tg,hS{l)) to indicate 
ignoring terms in hS{l) and in Tg = {symbol Psiy',v') which depends only on the 
s-derivatives of A(y', s), /i(y', s), p{y', s) up to order k}. 

Proposition 3.4. There is a bijective linear map W{y', s, rj') on the set o/3 x 3 
matrices which depends only on A, /i, p, but not on their normal derivatives, such that 

(3.19) X,{s){y\p')r^W{;s,-){DgX,+,{s)){y',Tj')) mod {Tg-^^+^\hS{l)), 
for any j < —1. 

Proof. Throughout the proof, we read x^ = s. First, we note that for j < — 1 
>^j{s){y',v') = -iD{x)Sj{x,r]'). 


From (3.15), we obtain 

SoS-i + S-iSo + D-^{R + R'^)S-i = -DgSo-D-^FoSo-D-^Fi mod {T°,hS{l)). 
Moreover, 

Fo = DgD mod (r°,h5(l)) 


and 


Fi = DgR^ mod (T°,hS{l)) . 


Hence, 

SoS-i + S-iSo + D-\R + R^)S-i = -DgSo - D-\DgD)So - D-\DgR^) 

= -D-^Dg{DSi+R^) 

= iD~^DgXi{s) mod (r°, h5(l)). 

Following the proof of Lemma 3.2, we find that A_i(s) satisfies 

(Q + pG)So^D-^X.iis) - X.iis)So = DgXois) mod (T^, hS{l)). 
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We note that 

Spec(5o) C C+, Spec(-(Q + pG)So^D-^) c C+. 

Hence, defining W{x,ri'){Y) as the solution X of 

(Q + pG)S^^D-^X - XSo = Y, 

we obtain 

,v') = W{-,s,-)DsXo{s){y',r]') mod (T°,hS{l)). 

For j < —2, Sj contains s-derivatives of A, /i, p up to order —j. Then, inductively, 
we get 

iQ + pG)So^D-^X,{s) - Xj{s)So = D,X,+iis) mod iT-^^+^\hSil)). 

Thus, we have proved the claim. □ 

We conclude this section by presenting the Riccati equation that A(s) satisfies: 
Corollary 3.5. Define 

(3.20) A(s) = i£)-iA^(s); 

A(s) satisfies, mod 0{h°°S{l)), the Riccati equation 

(3.21) /iiA,A(s) +Ji(s)A(s) + A(s)Ad(s)+A(s)2+F2(s) = 0 (0 < s < 1), 

where 

J^{s) = D-^{R + hFo), Kfis) =-D-'^R^, 

and 

F2(s) = -hDfiD-^ R^)-D-\R+R^+hFo)D-'^ R^+D-\hFi+Q+pG) + {D-^ R^f , 
with = s. 

Proof. This follows straightforwardly from (3.12) and (3.18). □ 


Invoking a forward Euler scheme to solve the Riccati equation, we obtain an 
approximate propagation of the boundary data into the interior of H, layer by layer. 
With the explicit reconstruction that will be presented in the next section, we obtain 
formally a layer-stripping algorithm for our inverse boundary value problem. The 
Riccati-type equation is expected to be highly unstable, especially for high frequency 
modes [18]. So the propagation of DN map will deteriorate. This reveals the ill- 
posedness of the problem of recovering the parameter in the interior. 

4. Reconstruction of the Lame parameters and density. 

We present the reconstruction of (A, p, p), as well as all their derivatives at = s 
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from the full symbol of A(s). We first consider the principal symbol of operator 
Lj fe 1=1 By the transformation rule of tensor, we have 

3 

N — (Y, 1 < *, fc < 3) = JNJ'^ 

(4.1) 

with N = (Y 1 < i, fc < 3). 

For any x near r(s), we choose a unit vector n{x) = (ni,n 2 ,n 3 ) € depending 
smoothly on x. Then any ^ = (^i, ^ 2 , Ca) ^ can be written as = qn{x) + m[x, C 
for some q G R and {mi,m 2 ,m 3 ) =■. m J- n. We define D = D{x), R = R{x,^), 
Q = Q{x,^) as follows, 


3 

D = iY, Cijkinjnv, I < i, k < 3), 

11=1 

3 

( 4 . 2 ) R={Y Cijkimjm; I < i, k < 3), 

jl=i 

3 

Q = (Y ^ijkimjmi; 1 <i, k <3). 

Jl=i 


We have (compare with (3.4)) 

(4.3) M = Dq^ + [R + {R)^)q -\- Q p. 

and a smooth factorization according to (3.5). More precisely, there exists a unique 
5'o = Sq{x,^) depending smoothly on a; G 17, and homogeneous of degree one with 
respect to ^ such that 

(4.4) M = {q- {So)*)D{q - So), Spec(^o) C C+. 

We let the direction of n{x) be aligned with the x^ axis. Using (4.1) we find that 

(4.5) M = {q- J{So)*J-^){JDJ^){q - {J^)-^SoJ^). 

Since the linear mapping defined by the matrix {j'^)~^ preserves the orthogonality 
m{x,^) _L n{x) and the length of n{x), we have q = 773 . We also have D = JDJ'^. 
Hence, by the uniqueness of factorization (3.5), 

(4.6) 5o = {J^)-^SoJ^. 

Combining this with (3.17) and the tensorial transformation R = JRJ^ of R, we 
obtain 

(4.7) Ao(s) = -iJ(U^o + if)J^- 

Due to the isotropy of the elasticity tensor, a rotation of coordinates {xi,X 2 ,X 3 ) 
does not affect the form D, Si, R and. Hence, for any x, we just assume ^ = (^', ^ 3 ), 
and n{x) = (0,0,1). Then m{x,^) = , so So{x,^) depends only on {x,^'). 
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and 


I) = diag(^,Ai,A + 2 Ai), A = , B = 


in which 

(4.8) 

with 


0 0 —ai 


a 0 0 


A = P 0 0 0 I F*, B = F 0 b 0 I F* 


-a2 0 0 


0 0 c 


Ciirr' C2icr' 0 
p = p( 0 =| 6 iri-' 0 

0 0 1 


(A + /x)|C'| 1 , MlfP + P 

Oil = ^ -, 02 = 70 : 1 , b=\ -, 

\/ fJ,{X + 2fi) 1 + 7 V ^ 


and 


.MlfP + P (A + Ai)2|^T 


1 + 7 


'(1 + 7)^ 


A + 2/i /^(A + 2^) 


a = 7 c 


7 = 


((A + 2 /i)|^'p + p){X + 2p) 


p{p\i'\'^ + p) 


We substitute = [I^'IjO]"^ and identify with |^'|; then, after some calculations, we 
find that 


Ao(s)(y,r?) 

= - iJ{b{x)SAx,0 + = JRoJ^ 

with 

Ai = (Af''\ 1 < i, fc < 3) 

ap 0 

0 bp 

ia 2 A/^(A + 2p) - iAIf'l 0 



ioi v^^(A + 2p) - i/r|^'| 
0 

c(A + 2/i) 


Since J is known and independent of A, p^ p, we only need to consider Ai for recovering 
A, p, p and their derivatives at = s. 

First step. We will recover A, p, p. Note that Aq^^^ {x, ^') = + pp. Then we 

can first get p and p as follows. Observe that 




-2 2 

= Cr, p , 


for any scaling constant cq > 0. Hence we set cq = 1 in the rest of this section. Then 
we find 


p = \J (x, \f2c^ ^)2 - Af {X, Cq ^ )2 
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and 

1 /; ( 22 ) / - 1\2 2 \ 

P ~('^0 (®> ^0 ) ~ P )■ 

P 

For A, first notice that we have 

^0 _ (A + 2/i + p)fj. 

(p + p)(^ + 2m)' 

Since we have already computed /i and p, we get 

A + 2/r + p _ 2 ^ P 
A F 2p A F 2p 

and then obtain A. 


Second step. We recover dX, dp,, dp of A, p, p, we first note that 

(4.9) 2pdp = d{X'^Q'^\x,V2co^f - 
from which we can recover dp. Then from 

(4.10) p5pFp5p = (9(Af^^(a;,Co^)^ 
we recover dp. Finally we recover dX from 


(4.11) 


dX + 2dp = d 


( ( M + 

^ Af)(FCo-i)2 



Final step. We recover higher order derivatives of X,p,p. Differentiating equations 
(4.9)-(4.11) k — 1 times, we obtain linear equations for d^p,d^X and d^p. The coef¬ 
ficients for them are the same as those for dX, dp, dp in (4.9)-(4.11). Thus we can 
recover d^p,d^X,d^p, using d^A.o{s){y,p) for j = 1,2,--- ,k, and d^p,d^X,d^p for 
j = 1,2, • • • , fc — 1. So we can recover all the derivatives of A, p, p recursively. 

5. Further implications. 

In this section, we show how to get a decomposition into incoming/outgoing waves 
via the factorization (3.5). In this section, we take T = oo and allow r to take complex 
values. We assume r is in the set 

IIo = {r G C; Kr := real part of r > 0}. 

For u G W((0, oo); D), we introduce the Laplace transform £ for t G IIo: 

poo 

{Lu)[x,t) = / e~'^*u{x,t)dt, 

Jo 

where {£u){-,t) G H^{£l). The inverse Laplace transform is given by 

2 I'-i+ico 

{C~^v){x,t) = — j e'"*^v{x,T)dT, 

ZTTl J -y_ioO 
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with 7 € IIq. 

Applying above Laplace transform to (1-1), we get 

3 

Mv = pr'^g^^Vk-h^ ^ V,= 0, 

j,k,l=l 

with h = 1 ^, f = 1 ^. Af can be viewed as a semiclassical pseudodifferetial operator 
with a small parameter h = 

We rewrite above equation up to the leading order terms in the following form 

(5.1) 

( hkv ) ~ ( -D-HQ + ^r^) -D-HR^R-) ) ( hL ) 

Denote M{x,T,ri) to be the principal symbol of Ai, as in (3.5), we have the 
factrorization 


(5.2) M{x,f,r]) = {r]3 - Sg {x,T,g'))D{x)(j ]3 - {x,T,r]')), 


for r G IIo. Here, similar to (4.8), 


(5.3) S+(x, f, r]') = 

and 


(5.4) So{x,T,r]') = JD^/^{A^-iB^)D-^/^J-\ 

where 


/ 0 0 -ai \ 

(5.5) A{x,T,g')=pi 0 0 0 P*, B{x,T,r]')=P 

\ -a2 0 0 / 


a 0 0 \ 

0 b O P* 
0 0 c/ 


with 


(A + At)|^'| 1 

\//r(A + 2fj,) 1 + 7 ’ 


a2 = 7ai, 


b = 


/mIC'P + 


c = 


1 

1 + 7 


'(1+7)^ 


A -l- 2^ 


(A + /r)^ieP 

/t(A + 2/i) 


a = 7 c 


and 


7 = 


/(A + 2/t)|fP+pr2)(A + 2/t) 


In all formulas for a, 6 , c, ai, a 2 , 7 , \/A is defined on C \ (—oo,0] with > 0. 
Indeed ai, 7 and b are well defined, and 5 R 7 > 0. Furthermore, we note that if Sr > 0, 
then } > 0 and $57 < 0, while 5(1 + 7 )^ < 0. Thus c is well defined for 

5r > 0. Similarly, we can verify c is well defined for 5r < 0. 
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We now show that 

(5.6) Spec(S'c;') C C+, Spec(S'o“) C C- 

for any r G Ho. The spectrum of Sq unified with the spectrum of Sq are the roots of 
det{M{x, T, rj)) as a polynomial in 773 for r G IIo. To obtain statement (5.6), we only 
need to show that there are no real roots of det(M(a;, f, 77 )) for any r G IIo. Then 
because (5.6) holds true for positive r, the eigenvalues of Sq or Sq cannot intersect 
the real line. Roots of det(M(a;, f, 77 )), which are same as the roots of det(M) in q 
satisfy 


+mI?T = 0 


or 


((A + 2/r)g^ + + pf^) (pg^ + (A + 2/i)|^'p + pr^) — (A + p)^g^|^^P — 0. 

For T G (0,oo), we have already concluded that g could not be real. For non-real r, 
if there is a real g satisfying the above equations, then 

= 0 , 


which is not possible. 

Let 

f iSQ-S+)-^SQ -iSQ-S+)-^ 

[-{Sq-S+)-^S+ {SQ-S+r^ 

and denote its inverse by P~^ which is given by 

s\)- 

Then we find that 

V -D-\Q + pGt^) -D-\R + R^) J \ 0 

It follows that 

satisfies 

{hDs - S+)v+ = 0 mod 0{hS{l)) 

and 

{hDs — Sq)v- = 0 mod 0{hS{l)). 

We can view u+ = C~^v+ {u- = C~"^V-) as representing incoming (outgo¬ 
ing) waves. This identification is justified by noticing that, for example, 7 ;+(-,r) = 
(£u+)(-,t) is exponentially decaying with increasing s for > 0 . 


0 


5n 


_ \P + \.O.t.. 
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We emphasize that the transformed DN map A, here, is different from the one 
introduced before, that is, 


=hCKTC-^- 

With the relation between the normal directive Dg and the DN map A, 
hDs = iD-^{K^ - R^) mod e)(/i5(l)). 


we have 


on boundary. We note that 

P = [So - s+)-^ 


V 

iD-i(A^ 


1 0 

0 -1 


P- 


mod 0{hS{l)) 


0 -1 

1 0 


with 


P-* 


1 5o- \ 

1 S+ 


where P~* = P~* for real r. 

After we have identified the elastic parameters on the boundary, we will then 
know R in (3.3). Then we can have a decomposition into incoming and outgoing wave 
constituents and u~ respectively) on the boundary using At and R. Seismic 
imaging (inverse scattering), array receiver functions, and tomography (also using 
free-surface multiple scattering) all rely on this decomposition. Discussion of seismic 
migration and inversion schemes based on this decomposition can be found in [6, 9]. 

Acknowledgement We thank the anonymous reviewers for very useful com¬ 
ments to improve this paper. 
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